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The Skyrme nuclear energy density functional theory (DFT) is used 
to model neutron-induced fission in actinides. This paper focuses on the 
numerical implementation of the theory. In particular, it reports recent 
advances in DFT code development on leadership class computers, and 
presents a detailed analysis of the numerical accuracy of DFT solvers for 
near-scission calculations. 

PACS numbers: 21.60.Jz, 21.10.-k, 21.30.Fe, 21.65.Mn, 24.75. +i 

1. Introduction 

In spite of successful applications in energy production and national se- 
curity, relatively little is known about the fundamental mechanisms of the 
nuclear fission process. Even today, the most widely-used theories of fis- 
sion rely on the macroscopic-microscopic approach to nuclear structure and 
semi-classical dynamics based on the Langevin equations. These methods 
are quite powerful but, in the long term, cannot be expected to yield the 
predictive power needed to understand the fission of very neutron-rich nuclei 
in the fission-recycling stage of the formation of elements, or to give enough 
accuracy for precise simulations of new generations of nuclear reactors. 

Already in the 1980ies, promising attempts were made to understand fis- 
sion in a microscopic framework based on the self-consistent nuclear mean- 
field theory with effective pseudo-potentials [H12]. At the time, the com- 
puting power was not sufficient for these approaches to compete with more 
empirical models, but these pioneer works yielded a lot of insight on the 
quantum mechanics of fission. In the recent years, the rapid development 
of leadership-class computers scaling to hundreds of thousands, and soon 
millions, of processing units has given us for the first time the computing 
power needed to successfully implement the full microscopic theory of fis- 
sion. In parallel, novel forms of scientific collaborations gathering nuclear 
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theorists, applied mathematicians and computer scientists have consider- 
ably improved our ability to utilize such large-scale machines [3]. 

The goal of this paper is to provide some tools to verify and validate 
DFT simulations of nuclear fission. In particular, we pay special attention 
to the problem of the numerical accuracy of DFT solvers at very large 
deformations. After a brief reminder on the theoretical framework, we recall 
some of the recent developments in DFT calculations on leadership class 
computers, then give a detailed analysis of truncation errors arising in DFT 
implementations using the one-center finite harmonic oscillator basis. 

2. Theoretical Framework 

The nucleus is described in the local density approximation at the Hartree- 
Fock Bogoliubov (HFB) approximation. The particle-hole channel is mod- 
eled by effective pseudo-potentials up to second-order derivatives in the 
density. In practice, this is equivalent to using zero-range effective inter- 
actions of the Skyrme type [H [5]. The results presented below were thus 
obtained with the SkM* parameterization of the Skyrme interaction [6j. 
The particle-particle channel is characterized by a density-dependent con- 
tact interaction with mixed volume and surface character [7]. An energy 
cut-off of E cnt = 60 MeV is used to reduce the number of quasi- particles in 
the definition of the densities. The HFB equations are solved in a one-center 
harmonic oscillator basis. 

In the DFT picture of nuclear fission, the HFB energy of the nucleus 
depends on an ensemble of collective variables q = (q±, . . . , qw). These can 
be, for example, variables describing the nuclear shape, excitation or spin. 
In this work, we considered as collective variables the expectation value (on 
the HFB ground-state) of the multipole moments Q\^- In practice, the axial 

Q20 and triaxial Q22, as well as the mass octupole Q30 and hexadecapole 
Q40 moments were considered. The collective space is thus four-dimensional. 
Expectation values of Qa^i will simply be denoted = (Qa^)- 

In the actinide region, the part of the potential energy surface relevant 
to nuclear fission spans a rather large range in deformations. The axial 
quadrupole moments runs typically from ~ 30 b in the ground-state to 
nearly 600 b at scission for symmetric fission; the octupole moment from 
to about 70 b 3 / 2 for very asymmetric fission (cluster radioactivity, see 
[8j [9] ) ; the hexadecapole moment from nearly ~ 3 b 2 near the ground-state 
to typically more than ~ 350 b 2 for symmetric fission. Assuming for sake of 
simplicity a uniform sampling of each degree of freedom, and a 1 b A//2 mesh 
size, the size of the collective space is more than 1.4 millions points, only 
for the axial collective variables. Adding triaxiality multiplies this estimate 
by another 2 orders of magnitude. 
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3. Large-Scale Potential Energy Surfaces 

In this section, we present the performance of our DFT solver, and 
provide a detailed analysis of convergence properties of our Skyrme HFB 
calculations in the case of 240 Pu. 



3.1. DFT Solvers on Leadership Class Computers 

All calculations were performed with the DFT solvers HFBTHO |10| 
and HFODD [11] . Both codes solve the HFB equations in the harmonic 
oscillator (HO) basis. The program HFBTHO assumes axial and time- 
reversal symmetry, while HFODD is fully symmetry-unrestricted. The two 
programs have been benchmarked against one another and agree within a 
few eV for an axial configuration |1U] . 
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Fig. 1. Performance of the DFT solver HFODD on the Titan supercomputer at 
the Oak Ridge Leadership Computing Facility in Oak Ridge. Each experiment 
measures the time of performing six full HFB iterations. The term 'deformation 
point' refers to a set of constraints on multipole moments, i.e. a point in the 
collective space. 



Owing to the block-diagonal structure of the HFB matrix induced by 
axial symmetry, the typical runtime of HFBTHO is of the order of the 
minute (depending on the type of nucleus, size of the basis and 'difficulty' 
of converging the HFB iterations). By contrast, it is of the order of several 
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hours for HFODD. In practice, HFBTHO is used as pre-conditioner for 
HFODD: for any given point q in the collective space, the HFB iterations 
are first solved with HFBTHO, and the densities at convergence are used to 
initialize HFODD. If the point of the collective space is axial, no additional 
iterations are therefore needed. 

Clearly, the very large size of the collective space together with the 
current runtime of the DFT solvers requires using today's most powerful 
supercomputers. A lot of effort was, therefore, devoted to porting our codes 
to leadership class computers, and ensuring that good scaling with the num- 
ber of processing units could be achieved. From a computational point of 
view, mapping the nuclear collective space in DFT is a naturally parallel 
problem. The code HFODD has, therefore, a hybrid MPI/OpenMP pro- 
gramming model, where points in the collective space are distributed across 
the MPI grid, and on-node multi-threading enables to take advantage of 
highly optimized linear algebra libraries. Figure [1] shows the performance 
of HFODD on the Titan supercomputer at the Oak Ridge Leadership Com- 
puting Facility. In this experiment, up to 300,000 processors were used in 
parallel. The slight degradation of the computing time is due to the orig- 
inal input/output backend, which has not been optimized and taxes the 
operating system at large scale. 



3.2. Numerical Accuracy 

Modeling nuclear fission requires to explore regions of the collective space 
with extreme deformations. The finite size of the HO basis may thus lead 
to a significant dependence of the results on the basis parameters. In our 
calculations, we only considered axially-deformed bases, characterized by a 
quadrupole deformation /3 and an oscillator frequency lvq. The maximum 
number of shells is denoted by iV max and the maximum number of states by 
-^states- In the case of a full spherical HO basis, the two are related through 
the well-known relation Abates = (AW + I)(N m3X + 2)(A 7 max + 3)/6. This 
relation is not valid anymore for a deformed basis. In practice, we introduce 
two cut-offs, one on the number of shells and another on the number of 
states. Note that the cut-off on the number of states is only the practical 
consequence of using a symmetry-unrestricted solver, for which the size of 
the matrices involved goes approximately as 2N^ ax . By contrast, the block 
structure induced by the built-in axial symmetry in HFBTHO would enable 
to consider all full shells up to N max . 

Basis dependence of the calculations thus comes from 4 parameters (i) 
the maximum number of shells N max , (ii) the maximum number of states 
Abates, (hi) the oscillator frequency uj and (iii) the basis deformation fa- 
in, principle, at every point q in the collective space, we should seek the 
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HFB solution that is the minimum in this 4-dimensional parameter space. 
Clearly, such a strategy is not sustainable even on today's largest comput- 
ers. Instead, one is bound to estimate truncation errors by exploring the 
parameter space locally, and extracting asymptotic expressions. 




Oscillator frequency 



Fig. 2. (color online) Convergence of the HFB energy as function of the oscillator 
frequency uq, for two configurations characterized by Q20 — 30 b (black squares) 
and Q20 — 300 b and Q40 = 120 b 2 (red circles). The deformation (3 is adjusted 
according to the formula (TJJ). 

As a first example, we show in figure [2] the dependence of the HFB 
energy on the spherical-equivalent frequency of the harmonic oscillator ojq, 
that is the frequency such that The deformation of the basis 

is fixed at each point according to the formula ([2]) . Two typical points in the 
collective space are considered, one near the ground-state with deformation 
Q20 = 30 b, one way past the second barrier on the descent to scission 
at Q20 = 300 b and Q40 = 120 b 2 . We note that the dependence on 
ojq is more marked at large deformations, and that the optimal frequency 
shifts toward smaller values as the deformation increases, which is consistent 
with the need to then include basis states with a larger spatial extension. 
Importantly, it is possible to extend this analysis and extract an empirical 
fit wo(<22o) giving the optimal basis frequency as function of the quadrupole 
moment of the collective point. In our tests, we found that the expression 

f 0.1 x Q 20 e-° m ^ + 6.5 MeV if |Q 20 | < 30 b m 
U ° \ 8.14 MeV if |Q 20 | > 30 b 1 ' 
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gives a reasonably accurate fit of the frequency as function of the quadrupolc 
moment. 




Basis deformation (3 

Fig. 3. (color online) Convergence of the HFB energy as function of the basis 
deformation (3, The oscillator frequency is adjusted according to the formula (TjQ). 
For the second configuration, an additional constraint on Q§ to Qqq — 150 b 3 was 
added. 

The dependence of the HFB energy on the deformation of the basis, for a 
fixed iV max and basis deformation, is illustrated in figure[3l Not surprisingly, 
the minimum is always obtained for basis deformation that are 'close' to the 
requested value of the axial quadrupole moment. Note that the dependence 
on deformation is rather marked. However, as for the oscillator frequency, 
it is a priori possible to obtain a fit /3(Q2o) such that the optimal basis 
deformation is chosen at point in the collective space. Our tests showed 
that the simple formula 

P = (2) 

provides a reasonable expression that remains applicable up to the largest 
values of Q20- 

Last but not least, we show in figured] the error induced by the trunca- 
tion of the number of states. Contrary to the previous two parameters of 
the HO basis, this truncation is imposed by memory limitations, and can 
not really be mitigated: for a given value of -/V max , the optimal number of 
states is always given by Abates = (AW + 1)(AW + 2) (AW + 3)/6, a 
number that can grow very large for large A^ax For example, at A^ax = 30, 
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Fig. 4. (color online) Convergence of the HFB energy as function of the number of 

States Abates ■ 



we have iV states = 5456. Taking into account the spin degree of freedom, the 
total number of basis states is more than 10,000, which implies that the size 
of the HFB matrix exceeds 20,000 x 20,000. At this time, it is not possible 
to handle in a reasonable time frame iterative processes involving dense, 
complex matrices in double precision of that size. As can be seen from fig- 
ure [U restricting A" states to manageable values around iV states « 1000-1200 
may easily lead to 2 to 3 MeV errors beyond the second fission barrier. 

To finish this section, we would like to emphasize two important conse- 
quences of basis truncation effects: 

• At constant truncation, the error increases with deformation, albeit 
not necessarily linearly. This is bound to have a very significant im- 
pact for, e.g., calculations of barrier penetrability, since the numerical 
precision is not the same at the entry and exit points, and errors of 
the order of the MeV can lead to orders of magnitude uncertainties 
on fission lifetimes; 

• Truncations magnify the impact of discontinuities in the potential 
energy landscape. In practice, calculations with very different basis 
characteristics initialized with similar density /wave- functions could 
converge to two very different points of the multi-dimensional PES. 
This effect is the reason why, in the lower panel of figure an addi- 
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tional constraint on Qqq had to be added: without it the calculation 
did not converge to the same point in the collective space at small and 
large basis deformations. 

4. Conclusions 

The nuclear energy density functional theory is currently the only viable 
option to achieve a microscopic description of nuclear fission. On-going de- 
velopment of leadership class computing facilities all over the world offer a 
unique opportunity to finally develop the nuclear DFT at very high preci- 
sion. In this work, we have discussed some of the numerical uncertainties 
associated with implementations of DFT in the one-center harmonic oscilla- 
tor basis. They clearly point to the need of developing bases that are better 
adapted to the extreme elongations characterizing the region near scission. 
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